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The atom-by-atom characterization of quantum gases requires the development of novel measurement tech- 
niques. One particularly promising new technique demonstrated in recent experiments uses strong fluorescent 
laser scattering from neutral atoms confined in a short-period optical lattice to measure the position of individual 
atoms in the sample. A crucial condition for the measurements is that atomic hopping between lattice sites must 
be strongly suppressed despite substantial photon recoil heating. This article models three-dimensional polariza- 
tion gradient cooling of atoms trapped within a far-detuned optical lattice. The atomic dynamics are simulated 
using a hybrid Monte Carlo and master equation analysis in order to predict the frequency of processes which 
give rise to degradation or loss of the fluorescent signal during measurements. It is shown, consistent with 
the experimental results, that there exists a wide parameter range in which the lifetime of strongly-fluorescing 
isolated lattice-trapped atoms is limited by background gas collisions rather than radiative processes. In these 
cases the total number of scattered photons can be as large as 10 8 per atom. The performance of the technique 
is related to relevant experimental parameters. 



I. INTRODUCTION 

The detailed characterization of non-linear many-body quantum dynamics is a major goal of ultracold atomic and condensed 
matter physics. In recent years there has been significant experimental progress in creating strongly interacting quantum systems 
within ultracold atomic ensembles (see the review article |lj| and references therein, and recent work including 

Hll). The 

strong interaction between neutral alkali atoms and laser radiation allows for the prospect that such strongly interacting quantum 
systems can be probed and investigated at the resolution of single atoms. This has been demonstrated very recently in two 
experiments 0|5|]; it may be anticipated that the excellent results from these experiments will stimulate further experimentation 
using similar techniques. 

The motivation to measure the spatial position of each atom in a strongly interacting quantum system is the desire to probe 
its structure and properties at an unprecedented level of detail. Strongly interacting quantum systems tend to exhibit much more 
complex dynamics than weakly interacting or linear quantum systems; this complexity leads to emergent behavior such as exotic 
phases of the quantum matter. Such behavior is often hard to simulate and fully understand on a theoretical level. Atom-by-atom 
measurements on these systems would allow much more detailed cross-checking of theoretical models against real data than is 
currently available, facilitating a much more complete understanding of these complicated dynamics. The substantial flexibility 
available in choosing the system Hamiltonian has led the idea that such systems may be used as 'quantum simulators' to probe 
various many -body Hamiltonians of interest in condensed matter physics J^H. 

In order to take images of the atomic distribution at single-atom resolution, recent experiments |0> HI 0] have used a deep 
far-detuned optical lattice to confine and mutually exclude the atoms as they are being measured. To detect the presence of 
an atom at a site with high confidence, the atoms must be localized to individual lattice sites throughout the duration of the 
measurement, despite substantial photon recoil heating. This means that the atoms need to be cooled as they are fluorescing; in 
current experiments polarization gradient cooling is used, with the cooling light collected to form the fluorescent signal. If atoms 
were to hop between lattice sites, not only would the spatial resolution of the signal be degraded, but there is a good chance 
that atoms would be lost from the measurement region entirely by undergoing a light-assisted collision with another atom, or 
J^j ' hopping to the edge of the lattice. Therefore, for efficient detection, hopping between lattice sites should be minimized while 
r\ , fluorescent scattering should be maximized. 

■ The purpose of this article is to model the physical process of polarization gradient cooling of isolated atoms in a deep far- 
detuned optical lattice in order to understand and support future experiments performed using this technique. In particular, the 
dependence of the atomic hopping rate between sites of the lattice is investigated as a function of experimental parameters. 

The available literature concerned with the polarization gradient cooling of confined atoms is limited to the study of atomic 
dynamics in dissipative optical lattices (for a comprehensive review, see flQID . and does not directly address the situation at hand. 
The polarization gradient cooling of ions has also been studied in separate work iTul [l2tl . 

In common with most theoretical treatments of laser cooling, this article uses a semi-classical Monte Carlo technique in order 
to model polarization gradient cooling in a computationally tractable manner. The current situation differs from previous work 
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in the context of dissipative optical lattices (e.g. |[l3l - [l5lP in the use of a separate far-detuned optical lattice to confine the atoms 
and a higher than usual atomic saturation parameter during the cooling process; these both help to maximize the fluorescent 
scattering to hopping rate ratio. 

This article also differs from previous articles by the use of a novel hybrid Monte Carlo-Master Equation (HMCME) technique 
to extract the quantity of interest — the hopping rate between wells — from the simulation. This is necessary due to the difficulty 
in extracting the hopping rate from the Monte Carlo simulation alone. The difficulty is due to the great disparity in time scales; 
while intra-atomic dynamics take place on a timescale of 10~ 6 s, it is desired that atomic hopping between lattice sites takes 
place on a timescale greater than 10 2 s. The direct Monte Carlo simulation of the hopping is therefore unfeasible. The approach 
taken uses Monte Carlo simulations of the short-time dynamics to construct a master equation rate model for the simulation of 
the long-time behavior. This is a novel approach for the simulation of ultracold atomic dynamics. 

The present article studies the dynamics of atoms undergoing polarization gradient cooling in a deep far-detuned optical 
lattice. While this is the basis for optical resolved-atom measurements of the spatial distribution of an atomic sample, other 
considerations are also important in optical resolved-atom experiments. Light-assisted collisions quickly eject atom pairs from 
multiply-occupied wells during the first small fraction of the measurement period; in current experiments the parities of the 
site occupation numbers are measured rather than the true spatial distribution of the atoms. Another consideration is that it is 
desirable that the design of the optical apparatus allows the resolution of each individual lattice site within the imaging plane. 
Furthermore, if the atomic sample is three-dimensional, scattering from atoms in the imaging plane should be distinguishable 
from the scattering from out-of-focus atoms. These considerations, while important, relate to the interpretation of the mea- 
surement data and to the design of the optical apparatus; as such, they may be divorced from the basic physical process under 
examination in this article, and are left to be discussed elsewhere lfl6ll . 

Looking beyond the two-dimensional measurements of the recent experiments JUBh, three-dimensional resolved-atom to- 
mographic measurements of atomic distributions are possible with suitable apparatus QSQ. Tomographic measurements require 
multiple exposures to probe the complete atomic distribution; this substantially increases the time needed for the measurement. 
This places correspondingly stricter requirements on the atomic fluorescence and hopping rates, and it becomes much more 
important to optimize the scattering rate from each atom without recoil heating leading to atom loss during the measurement 
period. 

The simulation developed in this article is used in a separate article by the author lfl6ll to inform the design of an experimental 
method which is capable of three-dimensional tomographic measurements of the position of each atom in strongly interacting 
quantum systems at half-wavelength spatial resolution. The article (Ref. Hill ) outlines a potential solution to the problems of 
resolvability and light-assisted collisions during fluorescent measurements of dense ultracold atomic systems. Furthermore, a 
solution to the out-of-plane scattering problem — the background noise generated by fluorescent scattering from atoms outside 
the tomographic section in a three-dimensional sample — is proposed, and is modeled using the simulation developed in the 
current article. 

In fact the fluorescent imaging of atoms undergoing in-lattice polarization gradient cooling is a measurement technique which 
has the potential to be used in a wider class of experiments than just the resolved-atom experiments described above. The ability 
to extract a very strong signal from each atom in a dilute sample may be useful in a variety of scenarios; for example, to make 
a time-of-flight (column-density) measurement on a very low atom-number sample which is unobservable using absorption 
imaging due to its dilution. Unlike for resolved-atom measurements on dense samples, very high spatial resolution is likely not 
needed, so simpler apparatus may be used. 

The structure of this article is as follows. A conceptual outline of polarization gradient cooling for an atom in an optical 
lattice is given in Section[n] with emphasis on how the lattice affects the cooling process. The system is described theoretically 
in Section [Hi] and the dynamical equations are rendered suitable for computation by the use of a semi-classical Monte Carlo 
wavefunction method. In Section [IV] the hybrid Monte Carlo-Master Equation (HMCME) method is developed in order to 
predict the atomic hopping rate between wells. In Section[V]the frictional force felt by atoms undergoing polarization gradient 
cooling in a far-detuned optical lattice is calculated and compared to the case of untrapped atoms. A calculation of the site 
hopping rate is presented for a specific set of parameters in Section [VT] The dependence of the site confinement time on the 
cooling configuration and system parameters is investigated in Section lVTIl and conclusions are drawn. 



II. POLARIZATION GRADIENT COOLING OF ATOMS TRAPPED IN A LATTICE POTENTIAL 

A. Polarization gradient cooling without the additional potential 

Polarization gradient cooling exploits the polarization gradients of light fields to cool atoms. We briefly review two spe- 
cific types of polarization gradient cooling which were analyzed in a seminal paper by Dalibard and Cohen-Tannoudji lfl7ll : 
the one-dimensional orthogonal linear (lin±lin) and orthogonal sigma (<x + cr~) configurations for two equal intensity counter- 
propagating laser beams. 
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(a) Without pinning lattice 




(b) In pinning lattice, minima coincide 
(0 = 0) 




(c) In pinning lattice, minima displaced 
(</> = n/4) 



FIG. 1: (Color online) Polarization gradient cooling in the ID lin ± lin configuration for an atom with J g — 1 /2. (a) Without 
pinning lattcice. The atom loses kinetic energy as it climbs the potential; near the top of the potential, it is optically pumped 
into the lower potential, losing energy over the complete cycle, (b) Confined by a pinning lattice (with spatial period Ai/2) with 
minima coinciding with the state-dependent lattice formed by the cooling light (<p = 0). The cooling is inefficient for atoms 
which have less than half of the lattice binding energy, (c) Confined by a pinning lattice (with spatial period Ai/2) with offset 
minima (<p = n/A). The low-energy curve crossing enables deeply-bound atoms to be cooled efficiently. 



For the lin JL lin configuration, the positive component of the electric field is 

E + oc cos (kz + <p) e~\ - i sin (kz + (f>) e+i (1) 

where e ± \ = + (1/ (x + iy). The light potential consists of alternating standing waves of <x + and cr~ light. For a ground-state 
atom with a spin of J g — 1 /2, the <x + standing wave couples more strongly to the mj — + 1 /2 state, and vice versa; the ground- 
state light shifts therefore oscillate alternately according to position (see Fig. [Tall. Furthermore, the cr + (cr~) light pumps atoms 
towards the ntj = +1/2 (mj = -1/2) state, which for negative laser frequency detuning, lies below the original state. The motion 
of the atom is therefore described by a loss in translational energy as the atom climbs the state-dependent lattice potential, 
followed by optical pumping from a peak to a trough in the lattice potential accompanied by a much smaller recoil momentum 
exchange (see Fig.[Ta|. This 'Sisyphus' process is dissipative. 

For the <x + cr~ configuration, the positive component of the electric field is 

E + oc cos {kz + <p) x - sin (kz + <p) y . (2) 

This polarization is always linear, but orientated at an angle which is proportional to z; it is convenient to think of a spatially 
varying atomic orientation axis parallel to the electric field, so the light is ji polarized everywhere. While the lights shifts are 
constant, they are not equal for ground state atoms with spin J g > 1 ; the Clebsch-Gordon coefficients are such that, with the 
laser tuned below resonance, the energy offset of the states increases monotonically with \ntj\. Optical pumping preferentially 
transfers population from the high \m j\ (higher energy) states to the low \m j\ (lower energy) states. As the atom moves, the new 
basis is different to that of the original basis. If the atom starts in a low \mj\ state in the original basis, the atom has a higher 
proportion of its density matrix in the higher \mj\ states in the new basis; these lie higher in energy, and the atom feels a net 
force opposing its motion. The optical pumping preferentially returns the atom into a low \ntj\ state at the new position, so the 
net force on an atomic trajectory is again dissipative. 
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B. One dimensional cooling in a lattice using orthogonal linear polarizations 

For an atom confined to a region smaller than a wavelength the nature of the polarization gradient cooling processes are 
modified. In the cases under consideration the atom is confined by an additional optical lattice potential, which shall be called 
the pinning lattice; the pinning lattice has an optical frequency far away from atomic resonances, so the mixing of the ground and 
excited states induced by the pinning lattice potential is negligible compared to the mixing induced by the cooling beams, and 
consequently scattering from the lattice potential is negligible. It is assumed that the pinning lattice potential is state-independent 
(or nearly so) for ground-state atoms, so the frequency detuning of the lattice is much greater than the width of the hyperfine 
structure of the excited state. 

In this section we consider a pinning lattice formed by counter-propagating laser beams; the lattice spatial period is A.i/2, 
where Ai is the wavelength of the lattice laser beams. The counter-propagating lin ± lin cooling light induces an additional state- 
dependent lattice potential with spatial period A c /2, where A c is the wavelength of the cooling light (in this article the cooling 
and lattice laser wavelengths will be similar but not equal, A/ m A c ). The pinning lattice potential is characterized by its depth 
and phase with respect to the cooling light (this phase will be a function of position if the lattice periods are not equal). It is 
assumed that the pinning lattice is substantially deeper than the potential induced by the cooling light. The pinning lattice may 
be in phase (Fig. [Tbl or out of phase (Fig. [TcJ with the lattice potential produced by the cooling light; these two situations have 
very different cooling characteristics. 

For the case of = 7r/4, there are equal intensities of <x + and cr~ light at the minima of the combined potential, with equal 
and opposite intensity gradients, so the potential minima for nij = ±1/2 are equal in energy, but offset a small distance from 
each other (Fig. [Tell. For laser light with a negative frequency detuning with respect to the atomic resonance, atoms are optically 
pumped from the upper potential to the lower potential. Therefore for <p = n/A atoms are polarization gradient cooled efficiently 
at low energies as the optical pumping direction reverses near the potential minima; Sisyphus-like cycles can be performed even 
by deeply bound atoms. 

For the case of <p — 0, the optical pumping changes sign half way up the pinning lattice potential (Fig. [Tbl so Sisyphus effects 
do not efficiently cool atoms with an energy below this point. Cooling can occur for atoms with energy greater than this point; 
however calculations will show (Figure [2c} that this is less efficient than for the <p = 7r/4 case. The equilibrium temperature of 
atoms undergoing lin J_ lin polarization gradient cooling with <f> — will be much higher than for those with (p-n/A. 

Another difference in the polarization gradient cooling of lattice-trapped atoms is that the additional pinning lattice potential 
rapidly changes the velocity and position of the atom. If the intensity of the cooling light is such that the optical pumping time is 
greater than half of a lattice oscillation period, the atom cannot undergo a full optical pumping cycle before the relative energies 
of the states are reversed by the motion of the atom (see Fig. [Tc]i, so the cooling efficiency will be decreased if the cooling 
beams are not intense enough, and the steady-state temperature of the atoms is expected to increase at sufficiently low cooling 
intensity (see Figs. [4al and [4bl. In comparison, consider that without the lattice the atom encounters similar light fields at a rate 
of approximately 2|v|//l c , while the lattice-bound atom encounters similar light fields at a rate of approximately ai osc /n. Thus the 
intensity below which cooling becomes degraded is affected by the frequency of the lattice oscillations; for deep enough lattices 
(i.e. fast enough oscillations) this intensity may be expected to be higher than is the case with no lattice present. 

A further point is that the polarization gradient cooling process relies on inducing a differential dipole potential on the sublevels 
in order to generate the cooling force. If the pinning lattice has a sublevel-dependent component then this modifies the cooling 
force, and can potentially disrupt the cooling process altogether. 

C. One dimensional cooling in a lattice using orthogonal circular polarizations 

The <x + cr~ configuration in one dimension, unlike the lin ± lin configuration, does not depend on the relative phase of the 
lattices as the light shift from the molasses is the same everywhere. A change in phase of the molasses is identical to a rotation 
about the axis of the laser beams; by a symmetry argument, it can be seen that the dynamics are unchanged in the direction 
parallel to the laser beams. (There is a redistribution of the components of spontaneously emitted radiation in the orthogonal 
directions upon such a phase change, but the effect of this on the cooling dynamics is minor). As with the lin _L lin configuration, 
the optical pumping should occur on a time scale shorter than half the period of oscillation for efficient cooling. 

D. Three-dimensional polarization gradient cooling of lattice-confined atoms 

Polarization gradient cooling can also take place if atoms are exposed to 3 sets of counter-propagating beams. The atoms 
are exposed to both polarization and intensity gradients in 3D molasses light; the polarization gradients cool the atoms by an 
admixture of the lin ± lin and cr + cr~ mechanisms [18]. 

Tightly confined atoms interact with the cooling light within a volume smaller than a cubic wavelength, so the net intensity 
and polarization that an atom experiences depends on all five of the relative phases of the molasses beams. As a consequence the 
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cooling efficiency along each of the three orthogonal directions at each pinning lattice site depends on the lattice site position 
together with these relative phases; while some atoms will be cooled efficiently, others will be cooled only weakly. Further- 
more atoms at different positions within the molasses light field will fluoresce at different rates; this can pose a problem when 
correlating fluorescent signal to atom numbers. 

This situation is unacceptable when the atoms are required, to a high probability, to be localized to each site for a long time. 
This article discusses two ways to overcome this problem. The first way would be to turn each set of counter-propagating one- 
dimensional cr + cr~ beams on sequentially, cooling the atoms along each direction in turn. The second way would be to introduce 
a small frequency difference between each pair of beams; the phases of the cooling light at any point, and so the molasses 
character, will then change continuously. This is the method used in recent experiments [:4!,i5j9j]. The two methods are analyzed 
in SectionslVllandlVlIl 

III. MONTE CARLO SIMULATION OF IN-LATTICE POLARIZATION GRADIENT COOLING 

This section outlines the Monte Carlo simulation of the polarization gradient cooling of atoms confined in an additional optical 
potential. 

The calculations presented in this article are performed for a single isolated atom, greatly simplifying the theoretical de- 
scription of the system. Reabsorption of scattered radiation is possible amongst an extended atomic sample; the influence of 
rescattering for a particular sample may be estimated by the ratio of the rescattered light to the incident light intensity at an atom. 
However, the atom number and density will be low in the situations under consideration, so this ratio will typically be much less 
than unity — the object of the measurement is to determine the position of each atom in a sample, and this will not be possible if 
there are too many atoms, or they are too closely packed. 

A. Simplification of the simulated system 

It is assumed that the additional confining optical potential has a frequency far away from the atomic resonance, so the mixing 
of the ground and excited states induced by the potential is negligible compared to the mixing induced by the cooling beams; 
consequently fluorescent scatter from the optical potential is negligible. As a further consequence mixing between internal 
states induced by the optical potential is ignored — of course, this mixing must be present to generate the potential, but in the 
far-detuned regime this mixing is much smaller than the mixing due to the cooling beams, and is ignored. 

In this article it is assumed that the additional optical potential is a three-dimensional optical lattice (the pinning lattice), 
with three pairs of beams frequency detuned from each other by many megahertz, so the effective potential is the sum of the 
one-dimensional potentials. The potential generated by the optical lattice is therefore 

V(r) = V (sin 2 kx + sin 2 ky + sin 2 kz) , (3) 

where Vo is a matrix dependent on the hyperfine structure of the transition; in general it is state-dependent, and includes non-zero 
off-diagonal ground-ground or excited-excited elements dependent on the optical polarization. 

The majority of alkali species have two ground hyperfine levels. Cooling takes place on a closed transition, i.e. \g, 1+ 1 /2) — > 
|e, / + 3/2), where / is the nuclear spin; however, the atoms also have a small probability of being excited to the \e,I+ 1/2) state, 
which can decay to the lower hyperfine state \g, I - 1/2). The atoms will eventually be optically pumped to this dark state. As 
usual, this behavior is to be prevented by the presence of 'repumping' light, i.e. light resonant with the \g, I - 1/2) — > \e, I + 1/2) 
transition. It is assumed that this repumping light is of sufficient power that the total residence time in the \g, I - 1/2} state is only 
a few times the natural lifetime of the atom in the exited states; consequently there is an extremely small population of atoms in 
the lower hyperfine state at any one time. This means that the time spent in the lower hyperfine level is much less than the period 
of vibration in the lattice, so negligible heating of the atom will occur due to differences in the optical potentials of the two ground 
hyperfine states. The repumping light is off-resonance by some gigahertz from the cooling transition \g, I + 1 /2) — > \e, I + 3 /2), 
and as such will only have a very small effect on the dynamics of the atoms within those states, which is ignored in the subsequent 
analysis. The population of the |e, / + 1 /2) excited state is also ignored. 

B. Unitary quantum dynamics 

The Hamiltonian describing the unitary evolution of the atomic state in the electric field of the cooling light in the rotating 
wave approximation is 

H = ~ \L V2 + V ° ( Si " 2 kX + Si " 2 ky + Si " 2 fa ) ~ HAP + ~T ( Es ( r)D s + E e( r ) D s) ■ ( 4 ) 

£=-1,0,1 
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The state vector contains 4/ + 6 components describing the \g, I + 1/2) and \e, I + 3/2) states. The matrix Vq is the state-dependent 
potential, and / the identity matrix. The matrices £)+ (D~) are the dimensionless raising (lowering) matrices with values given 
by the Clebsch-Gordon coefficients. These matrices are normalized so that 



£>:«; , (5) 

£=-1,0,1 



in which P is the excited state projection operator. The relative normalization of the local electric field vector E e {r) and the 
Raman frequency parameter Q. is arbitrary. The frequency detuning from resonance of the cooling light is A. 
The evolution of the density matrix p due to unitary processes obeys the von Neumann equation 

ih-p = [H,p]. (6) 



C. Relaxation mechanisms of the internal states 



The effects of the relaxation processes are determined from a trace over the 'environmental' degrees of freedom of the corre- 
sponding unitary processes. To account for recoil from spontaneous emission processes, recoil momentum terms are added to 
the standard result lfl9ll in a sum over all transitions and polarizations 



dt P 



= r 



-1(PP +P P) + 2 j 



COO 



(7) 



Here k is a unit vector centered on the origin, and the parameter f^ E ,( K ) describes the angular distribution of the spontaneous 
emission of photons with polarization <x for a particular combination of raising and lowering operators s and e'. The magnitude 
of the laser wavevector is denoted by Icr, and F is the angular frequency width of the excited state. The parameters /^{k) are 
derived in Appendix lAl 



D. Explicit retention of the excited states 



In many analyses of polarization gradient cooling (e.g. ITol [l7|p the excited states are adiabatically eliminated. This is 
appropriate for investigations that look to find the lowest temperature of atoms in optical molasses, as this usually occurs at 
low intensity; at low intensities, the population of the excited states is small, and adiabatic elimination is a good approximation. 
However, a high scattering rate is preferable for efficient signal extraction during the measurement process; furthermore, as 
discussed in Section HTBl the lowest temperatures of confined atoms tend to occur at higher intensity than for untrapped atoms. 
Therefore, the excited state populations are retained in the analysis in order to accurately simulate the situations of interest. 



E. The Wigner transformation 



In anticipation of taking the semi-classical approximation, the quantum dynamics of Equations[4]and[7]are expressed in terms 
of the Wigner function 

(8) 



W(r,p,f) = d 3 u(r+ ||p|r- |) e -^»/ fi . 
It is shown in Appendix |B1 that the complete dynamical equation in terms of the Wigner function is 

(~ + --v)w<r,p,t) = 

\ at m j 



2ni 

r 



■ (pw + wp) + r 



c- c' rr *J 



where 



d 2 KD;W(r,p - hk R K, t)D+f£(K) 



n 2 r 



Lm 



The Hamiltonian H is given by Equation|4] 



(9) 



(10) 
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F. The semi-classical approximation 

A full quantum treatment of the problem is computationally infeasible; instead, the semi-classical approximation will be 
used. In the semi-classical approximation the external degrees of freedom of the atom are treated classically, while the internal 
degrees of freedom are treated quantum mechanically. The semi-classical approximation is used in most theoretical treatments 
of polarization gradient cooling fiol [l3l [Til ITtI [l9ll . The approximation can be used in situations in which the coherence length 
of the atomic ensemble is much less than a wavelength; it relies on the dephasing influence of spontaneous emission to in effect 
'localize' the atoms. The accuracy of this approximation as used in the situation under consideration is discussed in Section 

ins 

Using the semi-classical approximation, the equations of motion are expanded in the small parameter e\ = hk/t/Ap, with 
higher orders in the expansion discarded Il9ll : here Ap is the momentum width of the Wigner function. For example, the 
expansion of the potential term on the right hand side of Equation|9]is 

V(r + s') = V(r) + s' ■ VV(r) + \- (V ■ V) 2 V(r) + ... . (11) 

This quantity predominately contributes to the integral in the region \s'\ < H/ Ap. The derivative of the potential V has approxi- 
mate magnitude V, so the series expansion is seen to be in the small parameter e\ ; it is terminated at appropriate order, which 
is chosen in this analysis to be the second order. 

This approximation leads to the equation (AppendixO 



in 



d 2 W d 2 V 



Zcrvr err \— < 

r m 2 ki ^ a 2 w . 

- - (PW + WP) + > TirttjD- „ „ Dt . 

e,e',i,j 1 1 J 

In this equation, the curly brackets are the anti-commutator, and the tensor rj EE >ij is given in Appendix lAl 



G. Conversion to Langevin form 



The semiclassical evolution equation as it stands (Equation \Y2\ still requires substantial computational power to simulate. 
Instead, the calculation will be restricted to a single trajectory; the distribution of atomic properties is then found by the sum 
over these trajectories. This approach is stochastic i.e. Monte Carlo in nature. The equations of motion for a specific trajectory 
in position and momentum space are found by substituting a semiclassical trial solution 



W(r, p, t) = w(t) 6°\r - f) 6 (i> (p - p) , 



(3), 



(13) 



into Equation fT2l This solution is valid in the limit in which both parameters e\ = Mu/Ap and ei = kRAp/mT are small, as 
discussed in Reference lfl9ll and in the previous section. The parameter ei ~ e{T /To, where To is the Doppler temperature, is 
much smaller than one in the situations considered in this article. 

Integration over external co-ordinates gives the equation of evolution for the internal co-ordinates 



at n 



-)-(Pw+wP) + Y J D- eP D + e 



(14) 



The equation of evolution of the external co-ordinates, e.g. f = (r), are found to be 

dn _ pi 
dt m 



dpi 
dt 



= fi = -TrAw 



dV(r) 



^{(n-n)(rj-rj)) = 



dt 

j t {(Pi-Pi)(Pj-Pj))=2D i j = tfklT—^Yj^ij^^K) ■ 



(15) 
(16) 
(17) 
(18) 
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The equation containing Djj describes the diffusion of the atom due to the atomic recoil. The diffusion term is a fluctuating 
Langevin force with zero mean, and is incorporated into the motion of a single trajectory in the Ito-Langevin equation 

dp i = f i dt + ^j2d~kdW k (19) 

in which dWj are independent, zero mean, Gaussian-distributed stochastic increments with variance dt. The quantities d^ are 
the components of the kth eigenvector of which is normalized according to its eigenvalue. 



H. Unraveling the Optical Bloch equations 



Equations[l5]to[T9]describe the motion of a single atom in the semi-classical approximation. However, the internal dynamics 
as described in[l4]are still ensemble-averaged, and are not appropriate to describe a single trajectory. The optical Bloch equation 
(Equation [Pit must therefore be unraveled into the stochastic evolution of a single wave function. The unraveling is chosen to 
be that of the quantum Monte Carlo wave function method (QMCW) lEoll . The underlying principle of the approach is to evolve 
the wave function using a non-Hermitian quasi-Hamiltonian with the addition of randomly-occurring discrete quantum jump 
processes. 

In the QMCW approach the following procedure is followed for each time step. Firstly, a random number < r < 1 chosen 
from a flat distribution and is compared to the quantity j = 1 - Y5t{ij/\P\4i). If r < j, 

m + 6t)) = \m) - istlv - l ^p) \m)- m 



n \ 2 

If the r > j, the atom undergoes a quantum jump (photon emission); photon polarization is chosen randomly according to the 
weights 

WDtDM) .... 

P £ = //ipi/\ ' (21) 

and the new wave function is given by 

W + 6t))=D;m)). (22) 

The new wave function requires normalization in either case. It can be shown that the ensemble of possible new wavefunctions 
satisfies the optical Bloch equation (EquationlT4l>. 



I. Computational methods 



The set of equations [T31 to l22l specify the stochastic evolution of a single semiclassical trajectory. An explicit third-order 
Runge-Kutta method is used to propagate the classical position and momentum of the atom. The coefficients for the method 
are chosen so that the intermediate evaluation times are h/3 and 2h/3, in which h is the time step for the external dynamical 
evolution. The internal components of the state vector are advanced at a constant time step h/3, and are propagated by the Cayley 
(split) form of the evolution operator. The new wave function is found efficiently by Gaussian elimination. 

The method outlined in this section was tested against simple analytic models and by comparison against previously published 
work looking at polarization gradient cooling in one-dimensional dissipative optical lattices Il3l - [l5ll . The results agree very well 
for intermediate saturation parameters. At high and low saturation parameters, the results differ somewhat; the retention of 
the excited states means the method presented here is more accurate for high saturation parameters (s>0.1), while the use 
of adiabatic elimination in the previous work enabled better statistics for lower saturation parameters (s<0.02). The Monte 
Carlo simulation was also tested against previous data for polarization gradient cooling of rubidium atoms in three-dimensional 
molasses ll2""ill ; for intermediate saturation parameters the predicted temperatures agreed to within 10%, which is around the error 
quoted for the previous work. 



IV. HYBRID MONTE CARLO-MASTER EQUATION ANALYSIS 



A. Comment on nature of the simulation 



The problem under consideration differs from previous analyses of polarization gradient cooling in a number of respects. 
The difference which poses the greatest computational challenge is the ratio of the time scales in the problem — between the 
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phenomena of interest, the jumping of atoms between wells, which we would like to take place over tens of seconds or greater; 
and the smallest time scale relevant to the problem, the period of the beat frequency between the laser and atomic resonance, 
which will be in the range of 10 ns to 100 ns. This means that, on average, one event of interest will occur every 10 9 to 10 12 time 
steps. Clearly a straightforward Monte Carlo simulation of these phenomena will be prohibitively slow. 

Adiabatic elimination is often used in simulations of polarization gradient cooling in order to increase the size of the smallest 
time step (and decrease the size of the state vector). However, this technique can only be used when the saturation parameter is 
much less than 1, which is not the case for the situations under consideration in this article (see Sect JIIIDb . In any case, the gain 
from this technique would be, at most, one or two orders of magnitude in the time scale ratio, i.e. interesting events would occur 
every 10 7 to 10 10 time steps; such calculations would still be very computationally intensive. 



In order to make the simulation tractable a method has been developed to extend the Monte Carlo simulation in order to study 
rare events. This technique bears similarity to previously developed techniques which are collectively known as the 'splitting' 



The basic idea is that a non-Markovian system can look Markovian when its dynamics are averaged over a sufficiently long 
time interval. An atom undergoing fluorescent scattering is clearly non-Markovian at timescales shorter than the scattering 
period; the deterministic state-dependent dynamics of the Bloch vector dominates, with stochastic dynamics (environmental 
coupling) playing only a minor role. However, when viewed over sufficiently long time intervals, it is known (from theoretical 
and experimental work) that quantities such as temperature can be assigned to an atom undergoing polarization gradient cooling; 
for such an assignment to make sense these quantities must be independent of the state of the atom at any one time, so implicitly 
it is assumed that the atomic dynamics, on long timescales, is described by stochastic Markovian dynamics. 

The hybrid Monte Carlo-Master Equation (HMCME) method presented in this section uses this idea to extend Monte Carlo 
simulation so that it is capable of predicting the frequency of rare events. The method is approximate; it uses the assumption 
that the system dynamics are Markovian when viewed over appropriately long time intervals. 



Firstly, a quantity should be found that is representative of the aspect of the system which is to be investigated; it should be 
a scalar time-dependent quantity which is a function of the state vector. It is chosen with two properties in mind; that it varies 
slowly and smoothly with time, and that the rare events of interest occur at values well separated from the usual values it takes. 
For confined atoms undergoing polarization gradient cooling, an energy-like quantity is appropriate; the atoms usually have 
energies well below that required to hop to a neighboring site. In the context of the subsequent analysis, the approximation will 
be used that the system performs Markovian dynamics — a random walk — along the energy axis. The Markovian approximation 
implies that the system is ergodic over long timescales. 

It is desired to label the system according to a discrete energy parameter, so a set of 'points' are defined at certain energies, 
{Ei , E%, . . .}. As the time evolution is simulated using a Monte Carlo method, the system will encounter these points repeatedly. 
(If the state vector of the system undergoes discontinuous jumps, the time and state vector of the system as it crosses a point may 
need to be extrapolated or approximated). The state of the system at any time will be classified by the last point encountered 
by the system. Each time the system encounters one of these points the state vector (the complete description of the system) is 
recorded; these records will be called 'start vectors'. 



A representative set of start vectors are required for each point, i.e. a set which fairly samples the true population of start 
vectors for that point. However, the run time of the Monte Carlo simulation should also be minimized. As the system encounters 
points at higher energies only very rarely, it would be prohibitive to take data for these points by means of direct Monte Carlo 
simulation. These two considerations need to be balanced. 

As a first stage, the simulation is run from some arbitrary starting point until it settles down into a steady state condition. From 
this time on the state vector is recorded every time the system encounters a point E,. A representative set of start vectors will be 
built up for points at low energy. 



B. An extended Monte Carlo analysis for rare events 




C. Classification according to energy 



D. Finding representative start vectors 
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E. Focusing the analysis on higher energies 



This coverage is now extended to points lying at higher energy. To do this a weight vv m - is assigned to each start vector; this 
is an estimate of how representative that start vector is of the general population of start vectors at all points. Initially an equal 
weight is assigned to each of the extant start vectors; the weights are normalized. The total weight of the start vectors at a 
particular point Wj is related, but not equal, to the probability of finding the steady state system in the environs of point ;. 

The nature of the simulation is now changed so that new start vectors ('daughter' start vectors) are to be found from existing 
start vectors ('mother' start vectors). To do this, a start vector is picked at random from all the possible start vectors according 
to its weight, and is evolved in time until a new point is reached. As it is assumed that the system is ergodic, the 'mother' start 
vector at the old point can be replaced by the new 'daughter' start vector at the new point (by assigning the mother start vector's 
original weight to the daughter and setting the mother's weight to zero), and the overall set of start vectors is still representative 
of the steady state dynamics of the system. 

Instead of the straight replacement of one mother start vector by one daughter start vector, the simulation is performed multiple 
times for each mother start vector; daughter vectors are found for the two neighboring points, together with the branching ratio. 
The mother vector is replaced by two daughter vectors, randomly chosen from those calculated, with one at each neighboring 
point; the weight of the mother vector is allocated to the daughters proportionately to the branching ratio. 

It is now possible to focus the simulation on higher energies in the confining potential where the atom is found more rarely. 
Once M mother vectors have been propagated for a point it is decided that enough data has been gathered to characterize the 
behavior of the system around this point, which is called 'full'. The start vectors at point are excluded from being propagated 
further using the Monte Carlo method; the choice of mother vector to propagation is now determined from the weights of the 
start vectors at all points which are not full. 

Although no further daughter vectors are generated from full points, weight is still added to full points from neighboring 
partially-filled points; if this goes uncompensated, it will lead to distortion of the distribution of weights between the points. 
Therefore as each mother vector is propagated and is replaced by daughter vectors, weight is reallocated amongst the full points, 
and from full points to partially-filled points, in order to preserve the relative allocation of weights between points. Details of 
this reallocation are given in AppendixiDl 

This method of extending the Monte Carlo simulation is not unique; in practice the method outlined was found to give the 
best balance between reproducibility, accuracy and computation time for the problem in hand. 



The method outlined in this section (and summarized in Appendix ID1 ensures that, as the simulation is run, it progresses 
from concentrating on the low energy region to the high energy region where the system is found only rarely. Along the way 
it has amassed data concerning the branching ratios between neighboring energy points (the probability that a system, starting 
at one point, will end up at either neighbor) together with the time taken to propagate from one point to the next. It is now 
simple to calculate transfer rates between the points. These transition rates are used in a master equation to calculate the relative 
populations of the atoms at various energies: 



where r mn = p m „/T m is the rate of transfer from m to n, t,„ the average time spent at m, and p mn the branching ratio. 

It remains to find the actual quantity of interest — the hopping rate between sites. The hopping rate of an atom at a particular 
point is directly extracted from the Monte Carlo data taken at that point; the total hopping rate is calculated as the average 
hopping rate for all points weighted by the relative populations at those points. The hop process is treated as a 'sink' for the 
atomic population, i.e. once the atom has departed the well it is removed from the simulation. 



A few subtleties were encountered when applying the algorithm to the polarization gradient cooling of confined atoms. Firstly, 
the potential energy of the system is ill-defined as the system has both quantum and classical properties. The usual energy 
measure in such a scenario is the expectation value of the energy Ek + (V(r)}, where r is the semi-classical position of the atom. 
However, this energy measure does not satisfy the two conditions outlined in Section lTVCI due to quantum jumps caused by the 
spontaneous emission of radiation, and by poor correlation between the value of the measure and the probability of site hopping 
(take the example of an excited state atom stationary at a lattice minimum). 



F. The master equation 




(23) 



III 



G. Implementation of the algorithm to the polarization gradient cooling of confined atoms 
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Instead, another energy-type measure was used which has the required properties. Using the kinetic energy as before, the 
potential energy was taken to be the lowest eigenvalue of the atomic system at the position of the atom. The time evolution of 
this quantity is continuous; furthermore, it correctly accounts for the combination of velocity and position needed to differentiate 
the main population of atoms from those which hop between sites. The choice is further justified by noting that atoms undergoing 
polarization gradient cooling in a sublevel-independent lattice tend to be optically pumped to this lowest energy ground-state 
sublevel (see Sect. Ill Bt . 

The energy spacing of the points was in practice determined by a compromise between precision and computation time; about 
25 points were used spread over the energy range to 1.3Vo, where Vo is the lattice depth. If more closely-spaced points were 
used, the computation proceeded faster, but at the expense of a larger spread in the prediction of the hopping rate between 
wells. This was attributed to insufficient diffusion of the state vector as it is propagated from one point to the next leading to 
observable 'clumping' of the simulated trajectories, and consequently undersampling of the phase space of start vectors at each 
point. Increasing the spacing between points, and so increasing the time the system spent at any particular point, decreased this 
trajectory clumping by increasing the influence of diffusive processes over shorter-time non-diffusive processes; consequently 
this decreased the variation in the predicted jump rate between runs of the HMCME simulation at the expense of increased 
computation time. Around 150 start vectors were propagated for most points, rising to ten times this number near the lip of the 
lattice potential. 

The accuracy of the algorithm presented in this section was tested by comparison with Monte Carlo simulations performed 
conventionally using the method described in Section[nI](which had been tested against previously published data — see Section 
1111 II I. There was excellent agreement between the HMCME method and the conventional Monte Carlo method throughout the 
region of overlap for all the data presented in this article (e.g. Figure |3bl. 



V. THE FRICTIONAL FORCE AND ENERGY LOSS RATE 

To investigate the effect of the pinning lattice on the polarization gradient cooling of ultracold atoms it is instructive to compare 
the form of the frictional force with and without the lattice present. Conventionally the frictional force is calculated for an atom 
on a constant- velocity trajectory lfl7ll ; however this does not appropriately describe the motion of atoms moving in an additional 
lattice potential. The average force and the average energy loss rate are instead calculated for an atom on a constant-energy 
trajectory. 

As has been discussed in Section lTV Gl there is ambiguity in the definition of an energy measure in the semi-classical situation. 
The measure used here for the potential energy (as in the rest of this article) is the position-dependent lowest-energy eigenvalue 
of the atom subject to both the pinning lattice and cooling light fields. The average force and energy loss rates are calculated by 
integrating Equations[14]to[T6]using a density matrix method. 



A. No additional lattice present 

The force and energy loss rates for polarization gradient cooling without an additional lattice present are presented in Figures 
[2a]and[2b] The dependence of the force on velocity for the <x + <x~ configuration is very similar to previous analyses |[l7ll . while 
the lin ± lin configuration looks somewhat different. 

The differences for the lin ± lin case are due to performing the calculation on a constant-energy rather than constant-velocity 
trajectory. For energies below 0.14 ftF the atom is confined to a single site of the optical lattice induced by the cooling light (see 
Fig- UK so the atom is not cooled efficiently. The atom experiences stronger cooling once it can move along the lattice. On the 
other hand, the induced dipole potential is spatially homogeneous in the cr + cr~ configuration, so atoms moving on trajectories 
which have constant energy also have constant velocity; results are found which are very similar to previous analyses. 

If an energy measure is used which does not contain potential energy contributions from the cooling light, the form of the 
lin ± lin velocity dependence alters to closely resemble that found in previous analyses (e.g. Ill7l0 . It is worth noting that in 
practice atoms will follow neither constant-velocity nor constant-energy trajectories. 



B. With additional lattice present 

The cooling force and energy loss rate for atoms in a pinning lattice of depth 1 AfiT are shown in Figures [2c] and [2d] The 
velocity plotted in Figure[2c]is the velocity of an atom at the minima of the lattice potential. 

The frictional energy loss profile for bound atoms is dependent on phase of the cooling beams at the pinning lattice site in the 
lin ± lin configuration but not in the cr + cr~ configuration. As discussed in Section HTTJ1 the cooling is very inefficient for atoms 
in a lattice with a <p — phase relative to the cooling light, but is efficient in the <p = n/4 case. 
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FIG. 2: (Color online) Polarization gradient cooling of 87 Rb in one dimension calculated on trajectories of constant energy, (a) 
Force versus velocity without the pinning lattice, (b) Energy loss rate versus energy without the pinning lattice, (c) Force versus 
velocity in the pinning lattice, (d) Energy loss rate versus energy in the pinning lattice. The vertical dashed red line indicates 

the depth of the lin _L lin cooling light lattice in Figure (b) and the depth of the bare pinning lattice (i.e. in the absence of 
cooling light) in Figure (d). The intensity of each molasses beam is 10.8 mWcirT 2 with frequency detuning -4F from the D2 
line. The pinning lattice potential with depth lAKT is generated by counter-propagating linearly polarized beams with 
frequency detuning +2000F from the Dj line. The parameter F is the angular frequency width of the D2 line. 



The character of the cooling is different for atoms which have an energy greater than the depth of the lattice. The cooling 
of unbound atoms is similar to the cooling of atoms without the lattice present — the lin ± lin configuration becomes phase 
independent, while the cr + <x~ configuration has greater cooling power. On the other hand, atoms with energies around the lattice 
binding energy experience reduced cooling power. This is due to time-averaging; these atoms travel slowly when near the lip 
of the potential, so they are cooled more weakly at that time, and the average force falls. These dips occur at slightly different 
positions due to the different depths of the overall potential in each case. 



VI. SIMULATED IN-LATTICE POLARIZATION GRADIENT COOLING 



A. Simulated scenarios 

As discussed in Section HTTJl three-dimensional molasses necessarily contains sub-wavelength scale intensity and polarization 
gradients. If an atom is confined so that it only samples the light field in a sub-wavelength region, the cooling efficiency is 
heavily dependent on the relative phases of the molasses light (as demonstrated in Section|VBj. Experiments requiring reliably 
efficient three-dimensional cooling for tightly confined atoms need to find a way to overcome this problem. 

Two experimental scenarios are simulated in this section; they differ in the method by which the phase problem is tackled. One 
method (which will be called the ID alternating configuration) uses a one-dimensional cr + <x~ cooling beam configuration which 
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TABLE I: Parameters for the scenarios discussed in Section[VT] The frequency detunings are given with respect to the closed 
F — 2 — > F' = 3 transition for the D2 cooling light and from the F — 2 — > F' — 2 transition for the D; pinning lattice light. The 

flash cycle time is defined as three times the cooling flash duration in any one direction. The frequency offset is the defined as 
the difference in frequency between the three sets of cooling beams (6 = v x - v y = v y - v z = (v A - v z )/2). The parameter T is the 
angular frequency width of the D2 line. The pinning lattice depth is given in Kelvin as E/Icb- 



Atomic species Species 87 Rb 

(794.98 nm) 
+2000r 
3.0 x 10 4 mWcrn" 2 
lAKT fjK) 
397 nm 

ID counter- propagating in each 
direction (intensities add) 
Linear 

D 2 (780.24 nm) F = 2 -» F' = 3 
-AT 

10.8 mWcirr 2 1.81 mWcirr 2 
Alternating ID <x + cr~ Offset 3D cr + o- 
18 jxs — 
18 kHz — 
— 6.1kHz 



TABLE EL Data from simulations using the parameters given in Table U The quoted site lifetimes are due to radiative heating 
only, and do not include the effect of background gas collisions. The quoted errors on the calculated parameters are the one 
standard deviation random errors and do not include systematic effects. The temperature calculations include kinetic and 

potential energy contributions. 



Cooling configuration 


Alternating ID a* a 


Offset 3D <r + o-- 


log 10 (Site lifetime/seconds) 


6.5±0.3 


6.6+1.6 


Mean site lifetime 


3.3 x 10 6 s 


3.8 x 10 6 s 


Mean temperature 


8.5±0.3 u.K 


10.3±2.0 u.K 


Mean scattering rate 


(1.89±0.001)xl0 6 s-» 


(1.41±0.11)x 10 6 s~ l 



is alternated between the three axes in turn. This can be realized in an experiment by using square-wave pulses to modulate the 
input to three acousto-optical modulators. In the other method (the 3D offset configuration), all six <x + cr~ cooling beams are 
used at once, with a small frequency difference between each set of counter-propagating beams. This sweeps the relative phases 
of the cooling light at each site. The initial values of the relative phases of the cooling light at the pinning lattice site were chosen 
to be random for each run of the HMCME simulation. 

The commonly used species 87 Rb is chosen for the analysis; similar results are expected to hold for other atomic species which 
are subject to polarization gradient cooling. The pinning lattice was chosen to be near the Di atomic transition to allow efficient 
filtration of the intense lattice light from the weaker fluorescence light during the measurement. The pinning lattice frequency 
detuning is chosen to be much greater than the hyperfine structure of the Di line, and the lattice polarizations chosen to be linear, 
in order minimize the potential energy differences between the ground-state sublevels (see Section lHBl and 181 12611 ). 



B. Discussion of the results 

The results of the HMCME simulation are shown in Figure [3] using the parameters given in Table Q] for the alternating ID 
cooling configuration. Figure [3a] gives the transition rates between points, while Figure [3b] gives the steady state populations at 
each energy point i.e. the proportion of the total population which last encountered that specific energy point rather than any 
other. 

The transfer rate to lower energy is greater than the transfer rate to higher energy for all but the lowest-lying points, i.e. on 
average the atom is cooled at all but the lowest energies. The difference in the rates increases with energy; this is to be expected 
as the average frictional energy loss increases with energy throughout this region (see Figure [2db. It is seen that the rate datum 
points have a small spread around the trend line due to the statistical uncertainty implicit in the Monte Carlo method. This leads 
to some variation in the site lifetime between simulation runs, and is reflected in the quoted statistical uncertainty on the lifetime 
(TableHB. 
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FIG. 3: (Color online) (a) Transition rates between different energy classes, as found from a single run of the HMCME 
simulation. Atoms are categorized according to energy according to the method described in Section|lV]and Appendix iDl 
Rates are given for transfer into the next higher or lower energy category together with the hop rate for transfer to a neighboring 

well. The data was generated using the HMCME method using the parameters of TableUin the alternating ID cooling 
configuration, (b) Histogram of population versus energy. The populations for the HMCME technique were calculated from the 
rates given in Figure (a); the Monte Carlo populations were found directly from a conventional Monte Carlo simulation. The 
depth of the potential, marked with the vertical black line, is 1.34SF, which is less than the bare pinning lattice (lAhF) due to 

the dressing with the cooling light. 



The 'hop' rate for a certain energy class is the rate at which atoms transfer from that energy class to a neighboring site; this 
dominates the behavior of atoms which have energy above the binding energy of the lattice. Only atoms with energy slightly 
above the lattice binding energy have an appreciable probability in this configuration of being captured by one particular site; 
atoms with larger energies will move between many sites before capture, or may leave the lattice. The effects of the jump process 
are seen in Figure[3b]in the sharp decrease of the population with energies above the lip of the lattice (it is assumed that, once the 
atom has hopped, it does not return to the original well). Note that the quantum mechanical tunneling of (deeply bound) atoms 
between lattice sites is negligible due to the substantial depth of the lattice. 

There is an excellent fit between the population predicted by the standard Monte Carlo simulation and that predicted by 
the HMCME simulation in the region in which the standard Monte Carlo simulation can be used. That the HMCME method 
simulates the system well in a known region gives good confidence that the HMCME method provides an accurate extrapolation 
of the dynamics to higher energies. 

Assuming ergodicity, the population data is fitted to a temperature once the density of states at the lattice site has been 
calculated using 

D ^ = 7T^3 f " V(r)d 3 r . (24) 

The fitted temperature population distribution (essentially the same when fitted to either the HMCME or Monte Carlo data) is 
plotted in Figure [3b] The fit is accurate in the low energy region, but the full HMCME method predicts more atoms in the tail 
of the distribution than would be expected at a given temperature. This is due to the description of the population in terms of a 
temperature being valid only when the frictional force divided by the diffusion coefficient is a linear function of velocity (see for 
example 12711 p66). The frictional force is in fact not linear — the slope of the <x + cr~ force decreases with increasing energy (Fig. 
[2cl ), so the population in the tail of the distribution is higher than what the temperature of the low-lying atoms would predict. 

The corresponding plots for the 3D offset configuration are similar to those displayed for the ID alternated cooling configu- 
ration. Calculated quantities — the temperature, fluorescent scattering rate and site lifetime due to radiative processes — for both 
configurations are given in Table HI1 The mean time before site hop is defined as the average time for an atom at the steady-state 
temperature at a pinning lattice site to leave that lattice site, with random phases of the cooling beams at the lattice site. Over 
a million photons are scattered per second in either scenario, with very long site lifetimes predicted. Of course, such long life- 
times will not be observed in experiments due to background gas collisions; the significance of these figures will be discussed in 
Section lyllll 

The quoted uncertainty in Table HH is the standard deviation on the mean of these quantities over a sample of 100 runs of the 
complete HMCME simulation; it is a combination of the uncertainty due to the random nature of the Monte Carlo simulation 
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together with a contribution which depends on the relative phases of the cooling beams at the pinning lattice site, as will be 
discussed in Section [VII A| The averages and standard deviations of the hopping rate are calculated in the logarithm throughout 
this article; this has a similar functional dependence to an average over the temperature. While this uncertainty on the mean 
hopping rate is large in percentage terms, it is negligible compared to the rate of background gas collisions. The greatest 
systematic uncertainty in the calculation is likely to arise from the use of the semi-classical approximation, which although 
necessary to make the problem tractable, is only approximately fulfilled; the parameter e\ = Hkf>/Ap, assumed to be much less 
than unity in this analysis (Sect. HITFT i. is 0.19 for the current example, with the atom scattering around 5 photons per oscillation 
period. It is worth noting that the high ratio of the dephasing event rate (i.e. the photon scatter rate) to the oscillation frequency 
makes a simulation based on transitions between pure quantum vibrational states problematic. To improve on the semi-classical 
approximation a fully spatially-dependent density matrix description is needed, which would be very difficult computationally. 



VII. LIFETIME BEFORE HOP VERSUS EXPERIMENTAL PARAMETERS 



The dimensionality of the phase space of the experimental parameters is fairly large, so a complete characterization of the 
system is not feasible within the scope of this article. Instead, the dependence of the hopping rate on individual parameters 
is investigated. One parameter at a time will be varied, with the other parameters of the system held fixed at the values given 
in Table U which have been chosen to give good performance. Both the ID alternating and 3D frequency offset molasses 
configurations, as discussed in Section [VI Al will be investigated. The relative phases of the cooling beams at the lattice site are 
not varied in this section (as this would greatly increase the number of computations needed). The effect of the relative phases 
will be discussed separately in Section [VII A| 

In fact, nearly all of the parameter sets presented in this section are predicted, in experiments, to have a site lifetime limited 
by the rate of background gas collisions rather than radiative processes; this will be discussed in Section [Villi The objective 
of this section, however, is a comparative study of the efficiency of polarization gradient cooling in a deep optical lattice; the 
quoted site hopping rate, although for most parameters unobservable in the physical system, in this case serves as a useful proxy 
measure for the cooling efficiency. 

The dependence of the temperature and the site hopping rate on the intensity of the molasses beams is shown in Figures [4a] and 
l4bl For both the ID alternating and 3D offset configurations the temperature passes through a minimum in the range investigated. 
At higher intensities, the temperature displays a linear dependence on intensity, the same dependence as for cooling without the 
lattice present lfl7i \WL The temperature increases markedly for low cooling intensities. It is suggested, as discussed in Section 
III Bl and in Reference that the behavior at low intensity is caused by less efficient following of the local equilibrium population 
by the atom in its motion when the optical pumping rate is lower, leading to a corresponding decrease in the efficiency of the 
cooling. The intensity at which the atoms are at minimum temperature is greater than for the case without the lattice present 
(for similar frequency detuning — c.f. for example OH). The time before the atom is lost from the site displays a complimentary 
relationship with the temperature, as expected. 

Figures|4c]and|4d]display the results of varying the molasses frequency detuning at constant saturation of the atomic transition. 
The saturation was kept constant in order to maintain an approximately constant scattering rate (and so an approximately constant 
optical pumping rate) as the frequency detuning is varied. To achieve this the intensity of the cooling beams were varied along 
with the frequency detuning A in order to keep the saturation parameter Q?/2(A 2 + T 2 /4) constant. 

The standard theoretical model of polarization gradient cooling for a J = 1 transition lfl7ll indicates that the temperature varies 
proportional to (a 2 /A) + A with a « 3F, with the transition point a between the two regimes dependent on the Clebsch-Gordon 
coefficients of the transition. In fact this form of dependence matches the data in the confined case reasonably well, albeit for a 
somewhat higher transition point. Again, the site hopping lifetime displays a complimentary relationship with the temperature. 

As expected, the site hopping rate depends strongly on the pinning lattice depth (Figures|4e]and|4f]i. If the depth of the lattice 
did not affect the cooling dynamics, the site hopping rate would be expected to exhibit an exponential dependence on the lattice 
depth. The line on the logarithmic plot is in fact not straight but has a negative second derivative, indicating that the cooling 
is less effective at larger pinning lattice depths; this is also seen in the slight increase of temperature with lattice depth. It is 
suggested that two effects may contribute to this. Firstly the oscillation cycle time decreases with lattice depth; the cooling 
efficiency is degraded, as discussed above, when there are too few spontaneous scattering events per oscillation cycle. Secondly, 
the ground sublevel-dependent component of the lattice potential increases with the lattice depth; this differential in general 
disrupts the cooling mechanism (as discussed in Section [LTTJb . disproportionatly affecting high-energy atoms which travel far 
from the lattice minimum. Nevertheless, in the given parameter range these effects are not enough to generate a net increase in 
the site hopping rate as the lattice depth is increased. 

The results of changing the spatial period of the lattice by altering the angles between the lattice beams are shown in Figures 



4g and|4h] In the range investigated (lattice spatial periods of 400 nm to 1000 nm), the temperatures remain approximately 



constant, as does the site hopping rate for the case of the 3D frequency offset configuration. The site hopping rate in the ID 
alternating configuration increases slightly with lattice spatial period, indicating somewhat less efficient damping of higher- 
energy excursions at longer lattice spatial periods. 
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Flash cycle time [[is] Frequency offset [kHz] 

(i) ID alternating (j) 3D offset 

FIG. 4: (Color online) The variation of the temperature (lower red line) and site lifetime (upper black line) against experimental 
parameters for atoms undergoing polarization gradient cooling in a three-dimensional optical lattice. All of the parameters 
apart from the parameter being varied are given in Table U except for the cooling intensity in Figures (c) and (d), which was 
varied along with the cooling light frequency detuning in order to maintain a constant saturation parameter. The phases of the 
cooling beams relative to the pinning lattice site were the same for all data. Five runs of the HMCME simulation were used for 
each data point; the sample standard deviation is displayed. For comparison, the dotted horizontal line indicates the 
approximate lifetime of the atom at the site due to background gas collisions. 



There is modest variation in the temperature and site hopping rate when the cycle time of the alternating ID configuration 
(Figure [4]} and the frequency offset of the 3D configuration (Figure [4j]i are varied (these parameters are defined in the caption 
of Table J]). The cooling efficiency decreases if the flash time in the ID alternating configuration is too short; this is to be 
expected, as the initial transient optical pumping period after the beams switch direction, during which the atoms are not cooled 
efficiently, takes a greater proportion of the total time when the flash cycle time is reduced. A similar effect occurs in the 3D 
offset configuration; if the frequency offset between the sets of beams is too great, the change in the nature of the cooling light at 
the lattice site can become fast enough to adversely affect the cooling dynamics. On the other hand, if the frequency offset (in the 
3D configuration) is too small, the lattice site is exposed to the molasses light with the less efficient cooling phases for a longer 
time each cycle; this encourages excursions to higher energies, and results in higher temperatures and decreased confinement. 

The effect of the pinning lattice frequency detuning was also investigated at constant lattice depth. In the direction of higher 
frequency positive (blue) frequency detuning there were only minor changes in the temperatures and site lifetimes (see Table 
Hill for A/ = +8000F). Lower positive frequency detunings were not studied due to the presence of features associated with 
near-resonant F = 1 — > F' — 1,2 (/' = 1 /2) transitions, which are not included in the model. 

A pinning lattice with negative (red) frequency detuning was also studied; the cooling frequency detuning was altered to 
account for the change in the energy of an atom at the lattice minimum. Although one may expect less efficient cooling due to 
the atom spending more time in positions with a greater energy differential between the ground sublevels, this in fact is not the 
case for the lattice studied (see Table [TTTJ). This is because the three-dimensional lattice studied has the linear polarizations of 
each constituent one-dimensional lattice orientated orthogonally, and consequently the energy differential between the ground 
sublevels cancels at the lattice potential minima. As a result, the temperatures and site lifetimes were not substantially different 
to the blue detuned case. However, the use of a red detuned lattice has a number of undesirable effects, such a substantial increase 
in the rate of spontaneous scattering of the lattice light, which increases the rate of unwanted transitions to the lower hyperfine 
state. 

The types of lattices used in the experiments described in articles J3] and JH were also studied for the cooling parameters 
listed in Table U (apart for the cooling detuning, which was altered to compensate for the effects of the red-detuned lattice). The 
results of the simulations (see Table ITTTb are similar to those quoted for the lattice of Section lVTl 



A. Comparison of ID alternating and 3D offset cooling methods 

The results of the simulations predict that both the ID alternating and the 3D offset cooling configurations produce excellent 
localization of a strongly-scattering atom at a pinning lattice site. However, the 3D offset configuration has a greater spread in 
the calculated temperatures, site hopping rates and photons scattering rates than the ID alternating configuration (compare the 
sample standard deviations in Table HI]). 

This variability in the 3D offset configuration is primarily due to the role of the five relative phases of the cooling beams. The 
[0, +6] frequency detunings of the three sets of beams periodically vary the character of the cooling light at each pinning lattice 
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TABLE III: Temperatures and site lifetimes for different lattice frequency detunings and spatial periods. The data in this table 
were taken (as with Figure|4|i for constant phases of the cooling beams. The quoted uncertainties are the sample standard 

deviations. 



Configuration 


ID alternating 




3D offset 


Parameters 


ThxK] lo{ 


;io(Lifetime [s]) 


T [uK] 


logio(Lifetime [s]) 


A ( = +8000r 


8.3±0.3 


6.0±0.2 


8.3+0.2 


7.0+1.0 


A; = -2000r 


8.1±0.1 


5.6+0.2 


7.7+0.3 


8.3+0.7 


A, = +8200r, period = 680 nm ([4]) 


8.7±0.2 


6.0+0.2 


7.1+0.2 


7.7+1.1 


1064 nm counter- propagating ([5]) 


8.4+0.2 


6.0+0.4 


6.8+0.2 


6.7+1.4 



site, but even so only samples a small portion of this 5-dimensional phase space, as only a single linearly independent phase 
parameter being varied in a cycle of frequency S. This means that changing the value of the four residual independent phases of 
the cooling beams does change the overall character of the light field at the lattice site, and this still affects the efficiency of the 
cooling at the lattice site, albeit to a lesser extent than if no frequency offset was used. 

To illustrate the residual variation in temperature and site lifetime, the worst confining set of phases for the 3D offset configura- 
tion were further investigated from the data set investigated in Section[VT](which comprises 100 runs of the HMCME simulation 
at randomly chosen initial phases). The temperature and site lifetime for this set of phases were found to be 9.8±0.4 u,K and 
10 A (2.2±0.3) s, showing substantially less confinement than the mean. 

The photon scattering rate is also a function of these relative phases, as can be seen in the standard deviation of the scattering 
rate (Table HB. which is 8% of the mean value; this variation is not good when performing accurate fluorescent measurements. 
The variation in the cooling and photon scattering can be suppressed using a more complex set of frequency detunings between 
the cooling beams; however, there is still likely to be some residual variation between lattice sites. 

It is therefore concluded that for the situations that were simulated, while the ID alternating and 3D offset configurations give 
similar temperatures and site hopping rates, the ID alternating configuration exhibits less variability between different lattice 
sites and cooling laser phase arrangements, and hence is more suitable for use for fluorescent measurements. 



VIII. DISCUSSION 

It is seen that for a wide range of parameters the loss rate from a site of the pinning lattice due to radiative processes can be 
made to be 10~ 2 s or lower. The loss rate of atoms from the whole lattice due to radiative processes will occur at a lower rate: 
to be lost from the lattice the atom would either need to encounter another atom and suffer an inelastic light-assisted collision, 
or travel to the edge of a lattice without being recaptured. Even so, the transfer of atoms between wells of the pinning lattice is 
unwanted as it acts to blur the spatial atomic distribution measurement. 

Regardless of the radiation scattering processes, the lifetime of atoms in the lattice is limited by collisions with room- 
temperature background gas particles, which inevitably lead to the loss of at least one atom from the lattice. The lifetime 
due to background collisions varies by experiment and species, but typically is around 10 2 s. If the lifetime due to radiation 
scattering processes is greater than around 10 3 s, the loss rate due to such processes will be negligible compared to the loss rate 
due to background gas collisions, and it can be said that the polarization gradient cooling is optimal. 

In comparison, a typical image exposure time used in experiments is around 0.5 sto Is 0>HI[^,|^]. Nevertheless, it is in 
general desirable to have the site lifetime as long as possible, which in practice means that it is limited by background gas 
collisions. This is because the site lifetime has a direct bearing on the fidelity of the measurement: for a 1 s measurement, 
1% of atoms are lost during the measurement if they have a 100 s lifetime, but around 10% are lost for a 10 s lifetime; this is 
significant if hundreds of atoms are present in the image. Furthermore, and more perhaps importantly for future experiments, 
three-dimensional tomographic measurements have longer exposure times than two-dimensional measurements, as more images 
need to be taken; consequently it is more critical that the site lifetime is as long as possible. 

In fact the results show that the radiative site lifetime is several orders of magnitude above 10 3 s for a wide range of parameters. 
This may be important for experiments as it adds to the robustness of the measurement technique. Consider that an additional 
source of heating is present, for example due to the reabsorption of scattered radiation; the increased heating will tend to diminish 
the site lifetime exponentially. If the parameters of the cooling light are chosen so there is a 'safety margin' between the radiative 
site lifetime and the background gas collisional lifetime, the measurement technique can be made robust to the extra heating. 

Another very important consideration is that in an experiment it may be useful to use additional laser beams, on top of the 
ones discussed in this article, in order to aid measurements on the sample. An example of this concept is discussed in a separate 
article iflrjll . where it is proposed to use an additional laser beam to be able to distinguish between atoms undergoing fluorescent 
scatter at different depths in a three-dimensional sample; this would be of clear benefit for tomographic measurements. The 
additional probing beam adds substantial extra heating to the atoms, and it is beneficial to use near-optimal cooling parameters 
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in order to diminish atom losses during the measurement. 



IX. CONCLUSION 



The polarization gradient cooling of atoms trapped in a deep off-resonant optical lattice was simulated using a hybrid Monte 
Carlo-Master Equation technique. The calculations indicate that there is a wide parameter range in which the lifetime at a 
lattice site of a strongly-scattering single atom undergoing polarization gradient cooling is primarily limited by collisions with 
background gas. This is consistent with recent experiments IH |H |H . 

The calculations suggest that for near-optimal parameters the radiative loss rate from a site is in fact orders of magnitudes less 
than the loss rate due to background gas collisions. This could be important for experiments, as it suggests that the technique 
can be made robust in regards to additional heating mechanisms not included in the present analysis, for example heating from 
rescattered light or from additional probing beams. 

It has been shown that by the use of an in-lattice polarization gradient cooling fluorescent measurement technique it is possible 
to extract large photon numbers of up to 10 8 from each atom during fluorescent measurements of a dilute atomic sample. This 
strongly supports the use of this technique when probing the spatial distribution of low atom-number samples of ultracold atoms. 
The author would like thank Christopher Foot (Oxford, UK), members of the Quantum Processes and Metrology Group (NIST, 
Gaithersburg, USA), and the NIST internal referees for their comments, and to acknowledge funding from the EPSRC (UK), 
Christ Church College, Oxford (UK), The Lindemann Trust (UK) and the Atomic Physics Division of NIST, Gaithersburg 
(USA). 



Appendix A: The angular distribution of spontaneously emitted radiation 

The intensity of the spontaneously scattered radiation is ( ll28ll p. 173) 

I(r) <x (E-(r)E + (r)) . (Al) 

The electric field of the scattered light in the far field is determined from the source-field expression ( ll28ll p. 328) 

E + a (r)cze a (e a -D-) (A2) 

in which the unit vector e a is the polarization vector of the radiation mode labeled by a. 

Dividing the intensity into two perpendicular polarizations, conveniently 6 and <p in spherical polar coordinates, the following 
relationships are obtained: 

Ig(r) cc {(0(K) ■ D + ) (0(k) ■ IT)) (A3) 
/,(r) <* ((fa) ■ D + ) (fa) ■ D~)) (A4) 

in which mode a has been identified by the unit wave vector k and the polarization. The expectation value runs over the internal 
degrees of freedom of the atom. By breaking down the dipole operator into the component basis, the intensity is expressed in 
the form 



I,(r) = ^fj]f^(D + e D- e ,). (A5) 



The coefficients f"L are listed in TablellVlwith the normalization chosen so that the total scattered power is ficoT{P) . 
The tensor r\ ee >ij of Equation[T2]is related to the functions f£ E , by 

Vss'ij = J] I ' **X* ' xj)JZ, (A6) 

with jfj being the unit vector in the i direction. The tensor rj ee >jj is given in Table fVl 
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TABLE IV: The parameters f°. of Equation|7]in terms of the polar angle 9 and azimuthal angle <p. 
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TABLE V: The tensor f] E efij of EquationfTSl 
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Appendix B: The Wigner transformation 



The transformation of the kinetic energy operator term of the unitary dynamics (Equation|Ui is covered in standard treatments 
of the derivation of the Wigner-Moyal equation 1I29I l30ll . The transformation of the potential energy operator term of Equation 
|4]proceeds straightforwardly using the position space representation of the Wigner transformation (Equation [SJ. For the non- 
unitary dynamics of Equation|7] the transformation of the terms containing the constant matrix P is trivial. The remaining term 
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of the equation is 

r Z f ^ e ' i "' t ' r ' D « PD ° e ~ ikRK ' r & (K) ■ (B 1 } 

For this term it is easier to use the Wigner transformation as expressed in the momentum representation, which is found to be 



p-f)e ir " /n . (B2) 



2 

The application of this transform onto the Expression lBll involves expressions such as 

e- ik * Kr \p) = \p-hk R K), (B3) 

which are evaluated by identification of the momentum translation operator. The matrices Df are constant in spatial coordinates, 
so the Wigner transformed term is 



'2/ 

ss'cr 



r > I d 2 KD- e W(r, p - Hk R K, r)D+/£ (*) ■ (B4) 



Appendix C: The semiclassical approximation and conversion to Langevin form 

The semiclassical approximation is taken by the termination of the expansion (Equation[TT|i of the right hand side of Equation 
[9] in terms of small changes in momenta. Terms containing polynomials in s are expressed in terms of the derivative of the 
Wigner function with respect to momentum, for example: 

V p f ^(r + i\ P \r-i)e-^-LJ d 3 s (r + S -\p \r - 0<T'>^ . (CI) 

The recoil term of Equation|9]is directly expanded in terms of small momenta: 

W(r,p + q', t) = W(r,p, t) + q'- V p W(r,p, f) + I (q' ■ V,,f W(r,p, t) + ... . (C2) 

Integrals are performed over the unit sphere, and a sum taken over polarizations. Terms containing odd powers of kappa vanish, 
leaving a residual second order term 

rn 2 k 2 R ^ a 2 w . 

-f?.^^' (C3) 

£,£ ,l,J J 

The tensor 7] EE >ij defined in Equation lA6l 

The conversion of Equation [12] to Langevin form is achieved by substitution of the trial solution of EquationQ~3] To find the 
equation of motion for the internal coordinates (Equation[T4b the substitution is performed directly and the integral taken over 
the external co-ordinates. To find the equations of motion for the external co-ordinates (Equations [T31 to [T8ll the trial solution 
is substituted and the appropriate expectation value taken (by integration over both internal and external coordinates). A key 
relation concerns the derivatives of the delta distribution 

Using this relation the equations of motion are found for f = (r) etc. The non-zero expectation values for the second moments 
of momentum are interpreted as a Langevin force. 

Appendix D: Summary of Hybrid Monte Carlo-Master Equation method 

A measure is chosen which is representative of the dynamics of interest (see Section lTV CI ), and a series of points are chosen 
on this measure. The value of the measure varies in time as the system evolves. The system is classified at any time according 
to the last point encountered by the measure. When the system's classification changes (when the measure encounters a new 
point), the state vector of the system is recorded; this is called a 'start vector' for the new point. 
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1 . (a) One or more initial standard Monte Carlo simulations are carried out. After the system settles into a steady state, the 

state vector is recorded every time the system encounters a new point; these form the initial population of start vectors 
for the rest of the HMCME simulation. This continues until enough start vectors have been found to characterize a 
few of the points on the measure (those where the system is most commonly found). 

(b) An equal weight is assigned to each of the initial start vectors. 

2. (a) A single start vector (the mother) is chosen randomly from the set of all start vectors at partially filled points with a 

probability which is proportional to the start vector's weight w mt h- 

(b) The mother start vector is propagated multiple times using the Monte Carlo method until the next point is reached. 
The resulting state vectors (the daughters) are categorized by end point; a single randomly chosen daughter vector is 
retained for each end point. The branching probabilities pij for the mother start vector are calculated. 

(c) Weights (wdt)j = PijWmth are assigned to the two daughter vectors. If one of the daughter vectors is at a full point, the 
daughter's weight is added to the combined weight at that point. 

(d) For each full point (with weight Wj) a new weight is assigned according to 

(Wj) = Wj - (Dl) 

with Wn the total weight at all partially filled points. The weight from this point is distributed to the two neighboring 
points. If a neighboring point (with weight Wit) is also full, it is assigned the new weight 

{W k ) new = W k + P -^l. (D2) 
W N 

where pj k are the calculated mean branching ratios from j to k. If the neighboring point k is not full the start vectors 
at that point which originate from point j are assigned the weights 

/ \ , PjkWmthWj W d , Hk 

[Wdt j^k ) = w dt j^k + — ^ (D3) 

where the summation 2 Wdtj-^k runs over all the daughter vectors of the point j at point k. The reallocation of weights 
described in this step is carried out so weight is reallocated from each full point. 

(e) A weight w m ,i, — is assigned to the mother vector. 

(f) If the number of start vectors propagated reaches the maximum allowed M at a point i, the point is called 'full'; the 
weights w, „ of all start vectors at this point are added to give W,- = £„ w,-,,,, and the mean branching ratios py are 
found for all relevant channels. 

(g) Part 2 is repeated until data has been taken for all points of interest. 

3. The results of the Part 2 are processed to form a master equation as discussed in Section lTV F| 

The reasoning behind the reallocation of weights from the full points (Step 2(d)) is to prevent a false buildup of weights at 
the full points; in the absence of this reallocation, weight is transferred from partially filled to full points, but not the other way 
around. To compensate for not propagating vectors from these full points, the weights are reallocated (Equations ID 1 1 to [D3l in 
the same proportion that would occur if start vectors were propagated from the full points. These equations are derived by noting 
that, when each start vector is propagated, a fraction w„„/,/W/v of the weight of Wn (the total weight at partially filled points, i.e. 
the total weight from which the start vector is chosen) is being reallocated; therefore, to compensate for not propagating start 
vectors from the full points, an equivalent fraction of the weight at full points Wj x (wmth/Wfi), for all full points j, should be 
reallocated according to the known branching ratios. 
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